# TODO: Add comment
# 
# Author: broken
###############################################################################
#import library
library(base)
library(raster)
library(gstat)


#path
inputPath= "E:/FIMO/APOM/data/pm/2014/"
outputPath="E:/FIMO/APOM/data/pm/rs/2014/"

#pm file

fileArray = list.files(path = inputPath,pattern=".tif$",full.names = FALSE,recursive=TRUE)
totalFile=length(fileArray)
#dem,nigh light, landcover
demFile='E:/FIMO/APOM/data/pm/codata/dem14.tif'
nlFile='E:/FIMO/APOM/data/pm/codata/nl14.tif'
lcFile='E:/FIMO/APOM/data/pm/codata/lc14.tif'

demCorArray=c(1:totalFile)
nlCorArray=c(1:totalFile)
lcCorArray=c(1:totalFile)

sCor=c(1:totalFile)
sRMSE=c(1:totalFile)
sRE=c(1:totalFile)

oCor=c(1:totalFile)
oRMSE=c(1:totalFile)
oRE=c(1:totalFile)

uCor=c(1:totalFile)
uRMSE=c(1:totalFile)
uRE=c(1:totalFile)

for(i in 1:totalFile){
	
	fileName=fileArray[i]
	
	#PM values
	pmFile=paste(inputPath,fileName,sep="")
	pmRaster=raster(pmFile)
	pm=values(pmRaster)
	corxy=coordinates(pmRaster)
	x=corxy[,'x']
	y=corxy[,'y']
	
	
	#dem values
	demRaster=raster(demFile)
	dem=values(demRaster)
	
	#night light values
	nlRaster=raster(nlFile)
	nl=values(nlRaster)
	
	#landcover values	
	lcRaster=raster(lcFile)
	lc=values(lcRaster)
	#cell id
	totalCell=length(pmRaster)
	cell = c(1:totalCell)
	
	#data frame
	table=data.frame(cell,x,y,pm,dem,nl,lc)
	newTable=table
	testTable=subset(table,pm<0)
	trainTable=subset(table,pm>=0)
	
	#save data frame to csv
	outFile=paste(outputPath,substr(fileName,1,30),"_frame_full.csv",sep="")
	write.csv(table,file=outFile)
	outFile=paste(outputPath,substr(fileName,1,30),"_frame_train.csv",sep="")
	write.csv(trainTable,file=outFile)
	outFile=paste(outputPath,substr(fileName,1,30),"_frame_test.csv",sep="")
	write.csv(testTable,file=outFile)
	
	#caculate corelation between aot,dem,nl,lc
	corelation=cor(data.frame(trainTable$pm,trainTable$dem,trainTable$nl,trainTable$lc))
	outFile=paste(outputPath,substr(fileName,1,30),"_corelation.csv",sep="")
	write.csv(corelation,file=outFile)
	demCorArray[i]=corelation[1,2]
	nlCorArray[i]=corelation[1,3]
	lcCorArray[i]=corelation[1,4]
	
	#caculate variogram
	empiVario=variogram(pm~1,locations=~x+y,data=trainTable)
		
	#sph fit
	sphModel=vgm(psill=200,model="Sph",nugget=35,range=4)
	sphFit=fit.variogram(empiVario,sphModel)
	
	#plot sph variogram to jpeg
	outFile=paste(outputPath,substr(fileName,1,30),"_vario_sph.jpg",sep="")
	jpeg(file=outFile)
	plot(empiVario,model=sphFit)
	dev.off()
	
	#-----------------------SIMPLE KRIGING-----------------------------------
	
	beta_mean=mean(trainTable$pm)
	simple_result=krige(id="pm",formula=pm~1,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y,beta=beta_mean)
		
	#save simple kriging result to csv
	outFile=paste(outputPath,substr(fileName,1,30),"_sk_result.csv",sep="")
	write.csv(simple_result,file=outFile)
	
	#edit tiff
	simplePMRaster=pmRaster
	simplePMValue=simple_result[,3]
	simplePMRaster[1:totalCell]=simplePMValue
	
	#save sk result to tiff
	outFile=paste(outputPath,substr(fileName,1,30),"_sk.tif",sep="")
	writeRaster(simplePMRaster,filename=outFile,format="GTiff")
	
	#cross-validation
	simple_result_3=krige.cv(pm~1,trainTable,sphFit,locations=~x+y,nfold=3,beta=beta_mean)
	outFile=paste(outputPath,substr(fileName,1,30),"_sk_3f.csv",sep="")
	write.csv(simple_result_3,file=outFile)
	
	
	#simple statis
	simple_cor=cor(simple_result_3$var1.pred,simple_result_3$observed)
	simple_rmse=sqrt(sum((simple_result_3$residual)^2)/nrow(simple_result_3))
	simple_re=sum(abs(simple_result_3$residual)/simple_result_3$observed)/nrow(simple_result_3)
	
	sCor[i]=simple_cor
	sRMSE[i]=simple_rmse
	sRE[i]=simple_re
	
	print("simple completed")
	
	
	#-----------------------ODINARY KRIGING-----------------------------------

     odinary_result=krige(id="pm",formula=pm~1,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
  
 	 
	#save odinary kriging result to csv
	 outFile=paste(outputPath,substr(fileName,1,30),"_ok_result.csv",sep="")
	 write.csv(odinary_result,file=outFile)
	
	 
	 #edit tiff
	 odinaryPMRaster=pmRaster
	 odinaryPMValue=odinary_result[,3]
	 odinaryPMRaster[1:7812]=odinaryPMValue
 
 	#save ok result to tiff
     outFile=paste(outputPath,substr(fileName,1,30),"_ok.tif",sep="")
     writeRaster(odinaryPMRaster,filename=outFile,format="GTiff")
 
 	#cross-validation
     odinary_result_3=krige.cv(pm~1,trainTable,sphFit,locations=~x+y,nfold=3)
     outFile=paste(outputPath,substr(fileName,1,30),"_ok_3f.csv",sep="")
     write.csv(odinary_result_3,file=outFile)
    

	#Odinary statis
	 odinary_cor=cor(odinary_result_3$var1.pred,odinary_result_3$observed)
	 odinary_rmse=sqrt(sum((odinary_result_3$residual)^2)/nrow(odinary_result_3))
	 odinary_re=sum(abs(odinary_result_3$residual)/odinary_result_3$observed)/nrow(odinary_result_3)
	
	 oCor[i]=odinary_cor
	 oRMSE[i]=odinary_rmse
	 oRE[i]=odinary_re
	 
	 print("Odinary completed")
	 
	 
	 #-----------------------UNIVERSAL KRIGING-----------------------------------

     universal_result=krige(id="pm",formula=pm~x+y,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
 
 
	 
	#save universal kriging result to csv
	 outFile=paste(outputPath,substr(fileName,1,30),"_uk_result.csv",sep="")
	 write.csv(universal_result,file=outFile)
	 
	 #edit tiff
	 universalPMRaster=pmRaster
	 universalPMValue=universal_result[,3]
	 universalPMRaster[1:7812]=universalPMValue
	 
	#save uk result to tiff
     outFile=paste(outputPath,substr(fileName,1,30),"_uk.tif",sep="")
     writeRaster(universalPMRaster,filename=outFile,format="GTiff")
 
 
	 #cross-validation
     universal_result_3=krige.cv(pm~x+y,trainTable,sphFit,locations=~x+y,nfold=3)
     outFile=paste(outputPath,substr(fileName,1,30),"_uk_3f.csv",sep="")
     write.csv(universal_result_3,file=outFile)
   
	 
	#Universal statis
	 universal_cor=cor(universal_result_3$var1.pred,universal_result_3$observed)
	 universal_rmse=sqrt(sum((universal_result_3$residual)^2)/nrow(universal_result_3))
	 universal_re=sum(abs(universal_result_3$residual)/universal_result_3$observed)/nrow(universal_result_3)
	 
	 uCor[i]=universal_cor
	 uRMSE[i]=universal_rmse
	 uRE[i]=universal_re
	
	 print("Universal completed")
	 
	 
	 #plot all to jpeg
     pmOrigi=pmRaster
     pmOrigi[pmOrigi<0]=NA
     outFile=paste(outputPath,substr(fileName,1,30),"_all_compare.jpg",sep="")
     jpeg(file=outFile)
     par(mfrow=c(2,2))
     plot(pmOrigi,main="Original")
     plot(simplePMRaster,main="SK")
     plot(odinaryPMRaster,main="OK")
     plot(universalPMRaster,main="UK")
     dev.off()
 
	
	str=paste("Complete",fileArray[i],sep=" ")
	print(str)
}



#export statis result
result_dataframe=data.frame(fileArray,demCorArray,nlCorArray,lcCorArray,sCor,sRMSE,sRE,oCor,oRMSE,oRE,uCor,uRMSE,uRE)
outFile=paste(outputPath,"_result_2014.csv",sep="")
write.csv(result_dataframe,file=outFile)
print("Finish All")
